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ABSTRACT 



Context. The current stellar atmosphere programs still cannot match some fundamental observations of the brightest stars, and with 
new techniques, such as optical interferometry, providing new data for these stars, additional development of stellar atmosphere codes 
is required. 

Aims. To modify the open-source model atmosphere program Atlas to treat spherical geometry, creating a test-bed stellar atmosphere 
code for stars with extended atmospheres. 

Methods. The plane-parallel Atlas has been changed by introducing the necessary spherical modifications in the pressure structure, 
in the radiative transfer and in the temperature correction. 

Results. Several test models show that the spherical program matches the plane-parallel models in the high surface gravity regime, 
and matches spherical models computed by Phoenix and by MARCS in the low gravity case. 

Key words. Stars:atmospheres 



1. Introduction 



. Atlas (lKurucdll970[ Il979h is a well-documented, well-tested, 

■ robust, open-source computer program for computing static, 
plane-parallel stellar atmospheres in local thermodynamic equi- 
librium (LTE). Since Atlas came to maturity in the 1970s, stellar 
atmosphere codes have progressed in a number of directions to 
include i mportant additional phy sics. For example, the Phoenix 
program (Hauschildt et al. 1999) includes advances such as the 
massive use of statistical equilibrium (NLTE) in place of LTE, 

, spherically symmetric extension in place of plane-parallel ge- 
' ometry, the inclusion of the dust opacities needed for brown 
dwarf temperatures, and the ability to include blast wave veloc- 
ity fields. 

These advances, while obviously moving toward greater re- 
ality, have not eliminated some persistent problems. For ex- 

■ ample, a detailed study o f Arcturus using Phoenix models 
' jShort & Hauschildt! l2003b found that theii- spherical NLTE 

models increased the discrepancy between the observed and 
computed spectra l irradiance. A similar analysis of models with 
solar parameters dShort & Hausc hildt 2005) also concluded that 
important opacity and/or other physics is still missing. 

This evidence that the significant improvements contained in 
the state-of-the-art programs have not achieved closer agreement 
with the observations of these fundamental, bright stars indicates 
a need to continue exploring additional physics. This is a valu- 
able role that Atlas can play because it is open source and freely 
available from Kurucz's web page {http://kurucz.harvard.edu), 
where anyone can download the code and use it as the start- 
ing point for studying other possible improvements. An ex- 
ample o f the advanta ge provided by this starting point is the 
work of lKupkal (1 19961) who explore d convection as represented 
by the full spectrum of turbulence ( Canuto & MazziteUilll991l: 
[Canuto et al. 1996) as an a lternative to the stand ard mixing 
length theory. More recentlv. lSbordone et alj (12007) have devel- 
oped a GNU-Linux port of Atlas for use in a range of applica- 



tions. In that spirit, we have developed versions of the Atlas pro- 
gram that replace the assumption of plane-parallel geometry by 
spherical symmetry. These codes, comparable to the LTE, spher- 
ical version of Phoenix or of the spherical version of MARCS 
dGustafsson et alj[2008h . can serve as the basis for the study of 
additional physics needed to understand luminous stars that are 
cool enough that NLTE effects are not dominant. Such continued 
studies are certainly necessary given the revolutionary achieve- 
ments of optical interferometry in imaging the surfaces of stars, 
thus providing powerful new observational tests of models of 
stellar atmospheres. 

2. From Atlas9 and Atlas12 to Atlas_ODF and 
Atlas.OS 

In addition to including a wide range of continuous opacities. 
Atlas also incorporates the opacity of tens of millions of ionic, 
atomic and molecular lines. The original treatment of these lines 
was via pre-computed opacity distribution functions (ODE) in 
the Atlas9 program. Later, Atlas 12 was developed to use opac- 
ity sampling (OS) of the same extensive line lists, and Kurucz 
continues to expand and improve the opacities for these codes. 
There are small parts of the two Atlas codes that must be dif- 
ferent to handle the different line treatments, but the majority of 
the codes are not affected by the line treatment, and these can 
be identical. However, as inevitably happens, small differences 
between the two codes develop over time; these can be seen by 
differencing the two source files. Therefore, before beginning 
our conversion to spherical geometry, we rationalized the two 
versions of Atlas to be identical in every way except where they 
must be different for the line treatments. Where there are differ- 
ences not associated with the treatment of the line opacity, we 
have adopted the Atlas 12 routine, under the assumption that it 
is likely to be newer and undergoing more active development 
than the older Atlas9 code. We have also converted our codes 
to current Fortran 2003 standards. To distinguish our programs 
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Fig. 1. A comparison of the temperature structures for two so- 
lar atmospheric models. The top panel displays the temperatures 
as a function of logj,)(Pgas). The solid line is the model com- 
puted with Atlas_ODF and the dashed line is the solar model 
atmosphere asunodfnew.dat from the Kurucz web page. The 
bottom panel shows the percentage difference between the tem- 
peratures of the two models, also as a function of the gas pres- 
sure. The differences are entirely due to the conversion to Fortran 
2003 and to using consistent, modern values of fundamental pa- 
rameters throughout the Atlas_ODF code. The numerical meth- 
ods are the same for the two models. 



from the Kurucz originals, we use the name Atlas_ODF for our 
version of Atlas9, and Atlas_OS for our version of Atlas12. 

To test the accuracy of our conversions, we computed a 
model of the solar atmosphere, starting from the input model 
asunodfnew.dat dated 2 April 2008 in the Kurucz direc- 
tory /grids/gridpQQodfnew. We leai-ned later that this model 
was com puted b y Castelli as part of her collaboration with 
Kurucz (iCastelli & Kurucz 20041) . We used the opacity distri- 
bution function bdfpQQfbig2 .dat located Kurucz's directory 
/opacities/d£p8Qf, which uses the solar iron abundance = 
-4.51. After 20 iterations, the flux errors are all < 0.2% and the 
flux derivative errors are all < 2.8%; most errors of both kinds 
are smaller than our output resolution of 0.1%. Figure [T] com- 
pares the temperature and pressure structures of the solar model 
computed with Atlas_ODF and the starting Atlas9 model. 

Starting from the Atlas_ODF code, we removed those com- 
ponents that used the opacity distribution functions and replaced 
them with the components needed to do opacity sampling. At 
this point we introduced several changes that are not present 
in Atlas12. Depending on the effective temperature of the star, 
Atlas12 adjusts the index of the starting wavelength, variable 
nustart, to eliminate wavelengths were the flux is negligibly 
small. In anticipation of future applications of the code, we have 
removed this adjustment to always begin at the shortest wave- 
length, independent of the effective temperature. This change 
has been propagated back to the Atlas_ODF version. Second, 
Atlas 12 always computes 30,000 wavelengths with a wave- 



length spacing equal to a constant spectral resolving power of 
10"^. We have modified this to be able to specify a spectral resolv- 
ing power < 10"*, and have the number of wavelengths adjust au- 
tomatically. We did this to test the dependence of the computed 
model on the spectral resolving power. Third, in assembling the 
master file of lines to be sampled. Atlas 12 uses a sorting rou- 
tine from the computer's operating system, which is outside of 
the source code. We have modified the subroutine prelinop, to 
perform this sort within the Atlas source code, making it self- 
contained. 

To test Atlas_OS, we computed a model of the solar atmo- 
sphere, again starting from the input model asunodfnew.dat. 
We used Kurucz's files lowlines, hilines and bnltelinesS 
for the atomic and ionic lines, as well as the molecular files 
diatomic, tiolines and h2of ast. After eliminating lines that 
did not contribute at the temperatures of the solar atmosphere, 
the number of sampled lines used in the calculation was about 
two million. Figure |2] compares the temperature structure of the 
solar model computed with Atlas_OS and the model computed 
with Atlas_ODF. The differences between the models using the 
two methods of including line opacity are comparable to the dif- 
ferences between our Atlas_ODF and the original Atlas9 shown 
in Fig.[T] Therefore, the joint agreements displayed in Fig.[T]and 
Fig- E] show that our Atlas_OS code matches the Atlas9 struc- 
ture. As an additional test of the two line treatments, we have 
used our modification to the opacity-sampled code mentioned 
earlier to vary the spectral resolving power. The opacity-sampled 
model shown in Fig.|2]used the default spectral resolving power 
of 30,000. Our tests found that repeating the calculation with 
spectral resolving powers of 10,000 and 3,000 produced almost 
no change in the resulting model structure. 

As an additional test, Kurucz (private communication) pro- 
vided us with a new Atlas 12 (opacity-sampled) solar model that 
we compare with our Atlas.OS model in Fig. [3] It is apparent 
that the agreement is very good. 



3. From Atlas.ODF to SAtlas.ODF 

The spherically symmetric version of the code, SAtlas_ODF, 
was created from Atlas_ODF. Plane-parallel models are labeled 
with the parameters effective temperature, T^ff, and surface grav- 
ity, usually given as \ogg in cgs units. For spherical models these 
two parameters are degenerate because the same value of log g 
is produced by different combinations of the stellar mass and ra- 
dius. Therefore, we have elected to use the three fundamental 
physical parameters luminosity, L,, mass, M,, and radius, R,. 
These can be supplied in cgs units or as multiples or fractions of 
the solar values, L^/Lq, M^/Mq and R^/Rq, where the values 
of Lq, Mq and Rq are available throughout the source code via 
a Fortran 2003 module routine. 

The radius, of course, will vary with depth in the extended at- 
mosphere. Therefore we have chosen to define the radius where 
the radial Rosseland mean optical depth, tr, has the value of 2/3 
because a photon has a probability of about 50% of escaping the 
atmosphere from that depth. Other choices could be ma de, such 
as Tr = 1 dGustafsson et al.| j2008). or T500 = 1 (iHauschildt et al.l 
119991). but these differences are nearly negligible. 

With the three basic parameters defined, there are three mod- 
ifications to the code: pressure structure, radiative transfer and 
temperature correction. 
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Fig. 2. A comparison of the temperature structures for two solar 
atmospheric models. The top panel displays the temperatures as 
a function of logj,)(Pgas). The solid line is for the model com- 
puted with Atlas_OS and the dashed line is the model computed 
with Atlas_ODF. The bottom panel shows the percentage dif- 
ference between the temperatures of the two models, also as a 
function of the gas pressure. The differences are entirely due to 
the way in which line opacity is included. Both codes use Fortran 
2003, the same fundamental parameters and the same numerical 
methods. 

3. 1 . Pressure structure 

Atlas9 and Atlas 12 both solve for the pressure structure in two 
locations. After reading in the starting model, the pressure struc- 
ture is solved from the initial temperature structure as a func- 
tion of Rosseland mean optical depth, T(tr), in the subroutine 
ttaup. After the first iteration the total gas pressure is computed 
by integrating the simple equation 



where m is the mass column density defined as 
dm — —pdr. 



(1) 



(2) 



and g is the constant gravitational acceleration in the plane- 
parallel atmosphere. 

In a spherical atmosphere there are three potential depth vari- 
ables: mass column density, m, radius, r, a nd radial Rosseland 
mean optical depth, tr. As discussed by iMihalas & Hummed 
([l974), there is no clear preference for any of these variables. 
Therefore, we have elected to adopt the radial Rosseland mean 
optical depth as our primary variable, and then use T{Tp^) to solve 
for the pressure structure by modifying the subroutine ttaup. 
This is done in the initialization of the calculation and for each 
iteration. 

The modifications to the subroutine ttaup include allowing 
the surface gravity to vary with the radial distance from the cen- 
ter of the star. 



g(r) = G 



M(r) 



(3) 



Fig. 3. A comparison of the temperature structures for two so- 
lar atmospheric models. The top panel displays the temperature 
as a function of log[Q(Pgas). The solid line is the model com- 
puted with Atlas_OS, and the dashed line is an Atlas12 model 
provided by Kurucz (private communication). The bottom panel 
shows the percentage difference between the temperatures of the 
two models, also as a function of the gas pressure. 

Because the mass of the atmosphere is negligible compared to 
the mass of the star, it is an excellent approximation to set 
M(r) - M,, giving 



(4) 



If the starting model is a spherical model, there will be initial 
values for r. If the starting model is plane parallel, we solve for 
r from the defining differential equation 



dr 



1 



pkRir) 



-drp 



(5) 



in its logarithmic form to minimize numerical errors. Here p(r) 
is the mass density, and kj^(r) is the Rosseland mean opacity, 
the sum of absorption, k, and scattering, cr, per unit mass as 
a function of depth, both of which are available from the in- 
put model. This solution begins by assuming an initial value for 
the atmosphere's extension, r(l)/7?,, after which the solution is 
performed using the B ulirsch-Stoer me thod given in Numerical 
Recipes in Fortran 90 dPress et al.lll996h . The fifth-order Runga- 
Kutta method, also from Numerical Recipes in Fortran 90, was 
also investigated, and is built into the source code, but we found 
the results from the two methods to be identical. After finishing 
the solution for r, the atmospheric extension, r(l)/Rt, is checked 
against the initial assumption. If the extension differs by more 
than lO"*" from the starting assumption, the starting assumption 
is updated and the solution is iterated until the extension con- 
verges to < 10"^. 

Using g(r), the equation of spherical hydrostatic equilibrium 

is 



dPtc 
dTR 



8(r) 
kRirY 



(6) 
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The solution begins at the upper b oundary by assuming an ini - 
tial value of k^il). Then, following 'Mihal as & Hummeif(ll974l) . 
we assume all properties, except pressure and density, are con- 
stant above r(l), and that these variables decrease with a constant 
scale height. This leads to 

PUD = ^(1)^- (7) 

^r(I) 

The gas pressure, Pg(l), is derived from the total pressure by 
subtracting values for radiation pressure, Pr{l), and turbulent 
pressure, Piil), if these are known. The gas pressure and the tem- 
perature are then used to interpolate an updated value for k^il) 
from the input model. This procedure is iterated until the upper 
boundary pressure converges to < 10"^. 

With the upper boundary condition established, Eq. |6] is in- 
tegrated for Ptot, again using the Bulirsch-Stoer method. At each 
step the gas pressure is found as described above, and the gas 
pressure and temperature are used to interpolate the correspond- 
ing value for the Rosseland mean opacity. 

This method of solving the hydrostatic equilibrium is also 
applicable to the plane-parallel atmosphere with g{r) - g and 
without solving for the radius. To test our implementation, we 
have incorporated the modified version of subroutine ttaup, 
with both the Bulirsh-Stoer and the fifth-order Runga-Kutta rou- 
tines, into Atlas_ODF and Atlas_OS. The maximum differ- 
ence between the Bulirsch-Stoer (or the fifth-order Runga-Kutta) 
method and the original Hamming method was less than our out- 
put numerical resolution of 1 part in lO** at all but two of the 72 
depth points. Therefore, the percentage difference between the 
methods is zero except at these two depths, where the differ- 
ences are only 0.021% and 0.014%. It is clear that the pressure 
solution is being done correctly. 

3.2. Radiative Transfer 

Atlas9 and Atlas 12 solve the radiative transfer using the inte- 
gral equation method. The complication introduced by a geomet- 
rically extended, spherically symmetric atmosphere is that the 
angle between a ray of light and the radial direction varies with 
depth. Numerous methods are availa ble for solving t his problem, 
of which we have chosen to use the .Rvbickil([T97ll) reorganiza- 
tion of the lFeautrie^(ll964 method. 

Following the approach described by 'Mihalas! (Il978h . we 
solve the radiative transfer along a set of rays parallel to the cen- 
tral ray directed toward the distant observer, as shown in Fig.H] 
A subset of these rays intersect the "core" of the star, defined 
as the deepest radial optical depth, which we usually set to be 
tr = 100. We sample the surface of the core using 10 rays. 
We tried both equal steps of jj. - cos 6 covering the interval 
1.0 < < 0.1 in steps of AyU = 0.1, which is shown in Fig.|4] and 
a finer spacing toward the edge of the core by distributing the 
rays as = 1.0,0.85,0.7,0.55,0.4,0.25,0.2,0.15,0.1,0.05. We 
found the results to be almost identical, so we chose to use the 
equal /i spacing. For the core rays, the lower boundary condition 
of the radiative transfer is the diffusion approximation. From the 
core to the surface we follow the convention used by Kurucz in 
his plane-parallel models by having 72 depth points. For the cen- 
tral ray these depth points are distributed eight per decade with 
equal steps of the AlogTR - 0.125 from logTR = 2 to -6.875. 
For the off-center core rays the steps will be different, depending 
on the projection. 

The tangent rays are those that terminate at the radius per- 
pendicular to the central ray, and which are tangent to a particu- 
lar atmospheric shell at that point, as shown in Fig.|4] The radial 
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Fig. 4. The geometry used for the ^ Rvbickil (Il97ll) method of 
radiative transfer The rays are distributed over the core in equal 
steps of 0. 1 in fi. 



spacing between the shells is set by the central ray, and these 
spacings define the impact parameters of all the tangent rays. 
With this geometry, we calculated values of fi at the intersection 
of each ray toward the distant observer with each atmospheric 
depth, and from these we compute the integration weights at 
each point over the surface of each atmospheric shell at every 
depth. The lower boundary condition for the radiative transfer of 
the tangent rays is the assumption of symmetry at the perpendic- 
ular radius. At the surface of the atmosphere the rays toward the 
distant observer have /i values that depend on the steps described 
above. When we want to use these surface intensities, such as to 
predict the observable center-to-limb variation, we map the com- 
puted I(fi) onto a specified set of yu-values using a cubic spline 
interpolation. 

This solution with the iRvbickil ( 1197 lb organi zation uses ex- 
actly the same equations as the original Feautrieil d 19641) method. 
Therefore, in the plane-parallel case, in which both can be used, 
the results must be exactly the same. To test this, we created two 
alternative radiative transfer routines for the plane-parallel codes 
Atlas_ODF and Atlas_OS, one new routine having the original 
Feautrier organization and the other having the Rybicki organi- 
zation, and we select which of the radiative transfer routines we 
want by using an input instruction at run time, holding every- 
thing else the sa me. The result o f the comparison is shown in 
Fig. |5] Using the peautried (1 19641) organization gives exactly the 
same result (the output files have zero differences), as it should. 
The tiny temperature drop at the top of the atmosphere shown 
in Fig. |5]is entirely due to the precise form of implementing the 
surface boundary condition. We explored different implementa- 
tions (changing one statement) that were logically equivalent, 
and we found the result shown to be the closest match to the 
original Atlas integral equation solution. The other implemen- 
tations gave the same temperature at the top depth of the atmo- 
sphere, T{ti), but have a slower convergence to the temperature 
derived using the integral equation method. 
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Fig. 5. A comparison of the temperature structures for two solar 
atmospheric models, both us ing the Atlas_ODF code. The solid 
line used the 'Rybicki' ("197 1 ) method to compute the radiative 
transfer, and the dashed hne used the original Atlas9 integral 
equation method. 



< 1 mK. Therefore, we have chosen to modify the Avrett-Krook 
temperature correction routine in the original Atlas codes to in- 
clude spherically symmetric extension. 

We start with the time-in dependent equ ation of radiative 
transfer in spherical geometry dMihalasll 19781) . 



dliv) l-f?dl{v) 



= k(y)[S(v)-I(v)], 



dr r d/i 

where S (v) is the source function given by 
S(v, T) = o-(v)J(v) + [1 - cr(v)]B(v, T). 



(8) 



(9) 



To be consistent with the approach in Atlas, we express Eq.[8]in 
terms of the mass column density (Eq.|2]i to obtain 



dliv) 1 - di(y) 



dm 



dfi 



k(y)[S(v)-I(v)]. 



(10) 



In general, Eq. \TU\ does not conserve the luminosity with 
depth because the atmospheric temperature structure is wrong, 
but we assume that small perturbations to the current structure 
will produce a constant luminosity with depth. That is, we intro- 
duce the perturbations 



m — mo 4- 6m, 
r - kq + 6r, 
r = To + 6T, 



(11) 
(12) 

(13) 



A note about the relative run times of the same code using the 
three different radiative transfer routines: the time per iteration 
using the Feautrier method is about half the time of the orig- 
inal integral equation method, while with the Rybicki method 
the time per iteration is about ten times longer than the original 
integral equation method. 

The difference in the execution time of the Feautrier and the 
Rybicki methods, which use the same equations, is due to the 
sizes of the matrices that must be inverted. The Feautrier method 
computes the radiation field for all fi values at each atmospheric 
depth, where the ^ values are the double-Gauss angles found 
to be superior bv ISikesI ( Il951h . The computing time to invert 
the matrices scales as the cube of the number of fi values, M^. 
We use three angle points because our tests using four and eight 
angle points are insignificantly diff'erent {AT < 1 K) from the 
three-angle solution, while the computing time is certainly in- 
creased. 

The Rybicki method computes the radiation for each individ- 
ual ray over all atmospheric depths. The number of depth points 
ranges from just two for the tangent ray that penetrates just one 
atmospheric depth, up to 72 depth points for the rays that reach 
the core. The need to invert these larger matrices causes the exe- 
cution of the Rybicki method to be longer. 

3.3. Temperature Correction 

Atlas9 and Atlas 12 perfor m the temperature correction using 
the Avrett-Krook method ( Avrett & Rrooklll963t lAvrettill964n 
modified to include convection (Kurucz 'l97d). While the other 
spherical atmosphere programs (Phoenix and MARCS) use their 
own m ethods to perform temp erature corrections, it is our expe- 
rience (iTvcner & Lesterir2002h that the Avrett-Krook method is 
extremely robust, capable of achieving temperature convergence 



and 

/(v) = /o(v) + 6Iiv), 



(14) 



where the subscript refers to the current structure. We also use 
the subscript for the current extinction, ko(v), and the current 
source function, Sq{v). 

The extinction and the source function are expanded to first 
order as 



k(v) = Ko(v) H • dm. 



dm 



and 



S (v) = Sq{v) H ■ dm H ^:::r- ■ oT. 



dm 



dT 



(15) 



(16) 



To simplify the notation further, we represent dko(v) / dm = k'^^iv), 
dSa{v)ldm = S}iv), and dSa{v)ldT = Sq{v). 
Using Eq.[TT|we can write 



dm ^ d<5m 
dmo dmo 
or 

dm — dmo(l + 6m'), 
where 

d6m 



6m' 



dmo 



(17) 



(18) 



(19) 



Therefore, the derivative in the first term in Eq. [TOlbecomes 

dl(v) _ dliv) 
dm dmQ 



(1 +6m')-\ 



(20) 
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In the second term in Eg. [TOl there is r^^ = (ro + i5r) Because 
we assume that 6r ro, we can perform a binomial expansion 
of keeping only the first two terms, to get 



1 



ro + 6r 



'"0 



6m \ 

pn)l 



(21) 



by using Eq.|2] 

Using these, the spherical radiative transfer equation, includ- 
ing perturbations, becomes 



-pfi d[Io{v) + 51{v)] l-fi^ 



1 + 



6m\ d[Io(v) + 6I(v)] 



1 + 6m' dnio ro \ pro j dp 

[koiv) + k'^iv) ■ 6m\[So{v) + 5;,(v) ■ 6m + So{v) ■ 6T 
- Io(v) - 6I(v)]. 

Clearing the 1/(1 + 6m') term and expanding Eq. l22l 
terms with second-order perturbations, gives 



(22) 
ignoring 



-pp- 



dlojv) 
dmo 



d6I(v) 1 

-p^^— — + 



- p^ dIo(v) 1 - p^ d6I(v) 



1 



6m' + 



dmo ro 
6m \ dIo(v) 



ro \ pro) dp 
+ [koiv) ■ 6m' + fc„ ■ 6m][Soiv) ■ 
+ ko(v)[S'o(v) ■ 6m + So(v) ■ 6T 
Note that the left hand side of Eq. 

dio(v) i-p^dioiv) 

OOTo ro dp 
and the right hand side has 

^o(v)[5o(v)-/o(v)], 



dp ro dp 

= ^o(v)[5o(y)-/o(v)] 
/o(v)] 

■dim. 

3] contains 



(23) 



(24) 



(25) 



and these equal each other because they are just the two sides of 
Eq. [To] for the current structure. Canceling these out of Eq. |23] 
leaves the first-order perturbation of the spherical equation of 
radiative transfer 



d6I{v) 1 - p^ d6I(v) 1 - p 
-PP— + — + 



6m' + 



6m \ dIo(v) 



dmo ro dp ro 
[ko(y) ■ 6m' + k'^iv) ■ 6m][So(y) - /o(v)] 
+ ko(v)[S'o(v) ■ 6m + S'oiv) ■ 6T - 6I(v)] 



pro) dp 



(26) 



The first angular moment of the first-order perturbation equa- 
tion is obtained by multiplying Eq.|26]by p and integrating over 
all p, to get 



dmo ro 



SJ(v)] 



- [ko(v) ■ 6m' + k'oiv) ■ 6m]Hoiv) - ko(v)6H(v) 

1 / , 6m\ 

- - Urn' + — \[3Ko(v) - Jo(v)]. 
ro \ pro ) 



(27) 



Dividing Eq.|22]by koiv) and integrating over all frequencies, we 
obtain 



1 



^o(v) 
- 6m' 



36K{v) - 8J(v) d6K(v) 
- -P 



ro 

Hoiv)dv 
6m 



6m 



Jo 



dmo 
' k',(v) 



dv = 



— \6m' H 

'"o \ pro 



koiv) 
[3Koiv) - Jo(v)] 
koiv) 



Hoiv)dv- 



dv. 



6Hiv)dv 
(28) 



We now assume that the correct choice of 6m will make the left 
hand side of Eq.|28]go to zero. That is, we assume that the pertur- 
bations of the radiation field, 6K and 6 J, vanish when the correct 
atmospheric structure is obtained . Thes e assumptions are equiv- 
al ent to the assum ptions used in iAvrettI (IT964). equation 25, and 
in lKurucj (1 19701) . equation 7.5. This leaves the right hand side 
of Eq. |28]as a differential equation for 6m, 



ao6m' + bo6m -I- co = 0, 
where 



o-o - Ho + 



[3Koiy) - Joiv)] 
koiv) 



bo 



-f 

Jo 



1 

ro 

k'oiv) 

Y^Hoiv)dv + 
koiv) pr^ 



prn Jo 



dv, 

[3Koiv) - Joiv)] 
koiv) 



dv 



and 



Co 



Jo 



6Hiv)dv ^6H = 9{-Ho, 



with 
9{ 



iAnr)^ 



(29) 



(30) 



(31) 



(32) 



(33) 



being the radially dependent Eddington flux that we need to 
achieve. The general solution to Eq. |29lis 



6m — - exp 



_ fM^dJ r 'J^e-rm^-dm, (34) 
J aoim) J J flo(m) 



where in is an integration variable. 

The correction for the mass column density found above has 
assumed that all the energy is carried by radiation. If the atmo- 
spheric temperature is cool enough, significant amounts of en- 
ergy can also be carried by convection in the deeper, less trans- 
parent levels of the atmosphere. Atlas calculates the convective 
energy transport by the mixing length approximation. The equa- 
tions in Kuruczl ( ll970l) do not contain the radial variable explic- 
itly, but they do contain the surface gravity, g, which now varies 
with r. However, the implementation of those equations replaces 
g in terms of the total pressure, which now implicitly includes 
the geometry. Therefore, there is no need to modify the original 
Atlas code to include convection in the spherical temperature 
correction, and Eq. |29] remains the same, with the addition of 
convective terms in the coefficients ao, bo and co as follows: 



ao 



bo 



-- //o(rad) + 
-I- //o(conv) 

r ^//„ 

Jo ^o(v) " 



1 r 

ro Jo 



[3^:o(v) - Joiv)] 



3V 



2(V- 



(v)dv -H 



dT 1 

+ Hoicom) — - 
dm T 



1 - 





koiv) 


7 


(- 


Vad) 




1 




P^t 


Jo 




9D 



dv 



D 



£> + V - V,d 
[3Koiv) - Joiv)] 



(35) 



^o(v) 



dv 



-I- 



dV 



2(V - V,d) dm 



1 + 



D + V - V,d 
D 

D + V- V„a 



and 

co^-H- 



Hoirad) - //o(conv). 



(36) 



(37) 
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In Eq.[35]and Eq.[36]the V is 
dlnr 



dlnf 



(38) 



D is from lKuruczl (1 1 91(f), and in Eq.[36]the Hq{v) is the radiative 
flux. Solving Eq. [34] for Sm, using the coefficients in Eq. [35] 
Eq.[36]and Eq.[37] the corresponding temperature change based 
on conserving the flux is 



dT 
am 



(39) 



Atlas uses two additional temperature corrections near the 
surface, where the flux error loses sensitivity. One correction is 
based on the flux derivative. Because this correction applies high 
in the atmosphere where the gas is quite transparent, radiation 
will carry almost all the energy, and it is a good approximation to 
ignore convective energy transport. The zeroth angular moment 
of the spherical radiative transfer equation (Eq.fTOli is 



dm 



^r'k(v){J{v)-S{v)l 



(40) 



Replacing J{v) by A[5(v)], expanding the Planck function in 
Siv) in terms of T, integrating over frequency and retaining just 
the diagonal terms of the A operator, the resulting temperature 
correction becomes 



\_ f d(r^'H) _ d[r-H(rad)] \ 
r- { dm dm J 



ir^wiT^m^^[i-(v)] 



dB(v,T) 
dT 



dv 



(41) 



The terrnAdiajs approximated by the plane-parallel expression 
given in lKurucj ( [l97 0), assuming it has minimal dependence on 
the geometry. 

A third temperature correction is used in the original Atlas 
code to smooth the region of overlap between the first two cor- 
rections. This is 



SUlf 



T6H 
477' 



(42) 



and this is retained here. The total temperature correction is 
therefore. 



6T = 6Tnux + STf, + ST, 



surf ■ 



(43) 



4. Comparisons between SAtlas_ODF and 
Atlas.ODF 

One test of the validity of the spherical code is to compute 
a spherical solar atmosphere, which should be nearly identi- 
cal to the plane-parallel model. For both models we used the 
Kurucz file asunodfnew . dat for the starting model and the 
file bdfp00fbig2 . dat for the opacity distribution function. To 
eliminate other possible sources of diff'erences, we used the 
Bulirsch-Stoer solution to solve for the pressure structures and 
the Rybicki method for the radiative transfer in the plane-parallel 
model as well as for the spherical calculations. The spherical 
model used the atmospheric parameters Lq - 3.8458 x lO-'^' 
ergs/s, Mq = 1.9891 x 10^^ g, and Rq = 6.95508 x 10'° cm. 
These correspond to T^s = 5779.5 K and log^ = 4.43845, 
which are slightly diff'erent from the canonical values used by 
Kurucz. Therefore, we computed the plane-parallel model with 
the consistent values of T^^ and log g. The comparison is shown 
in Fig. [6] 




Fig. 6. A comparison of the temperature structures for two solar 
atmospheric models. The top panel displays the temperatures as 
a function of log[o(Pgas). The solid line is computed using spher- 
ical geometry and the dashed line is using the traditional plane- 
parallel geometry. The radiative transfer and the pressure struc- 
tures of both models were computed using the same numerical 
techniques. The bottom panel shows the percentage difference 
between the temperatures of the two models, also as a function 
of the gas pressure. 



The |Ar| is < 0.25% until log[oPgas < 2, where the tem- 
perature of the spherical model begins to trend lower than the 
plane-parallel model. The dip in the temperature difference down 
to -2.4% is due to the kink in the temperature structure of the 
plane-parallel model. This feature was discusse d earher i n con- 
nection with Fig. [5] However, in this case the iRvbickil (il97ll) 
method is used to compute the radiative transfer in both mod- 
els, using the same surface boundary condition in both codes. 
Therefore, this difference cannot be due to a coding difference 
between the two routines. This feature might be due to the num- 
ber of rays used in the calculation of the radiative transfer. The 
Rybicki solution for the plane-parallel code uses three rays for 
each depth, whereas the same method in the spherical code uses 
a; 80 rays for the layers approaching the surface. Perhaps this 
finer griding produces a smoother temperature profile in these 
layers. There is, however, no physical significance to the tem- 
perature differences near the surface because these lay ers are 
located in the solar chromosphere ( Fontenla et al J 12006). well 
above the temperature minimum, where other physics is com- 
pletely dominant. 

A test where larger differences are expected is for the 
coolest model (Teff = 3500 K, log^ = 0.0) in the grid 
/grids/gridp00od£new/ap00k2od£new. dat (computed by 
Castelli) on the Kurucz web site. Because T^^ and \ogg re- 
ally represent the three parameters L, M and R, there is a de- 
generacy that must be broken. To do this, we have assumed 
the star has M - I Mq, which leads to L» = 3690 Lq and 
R = 166 Rq. The comparison of the atmospheric structures 
is shown in Fig. [7] While the plane-parallel model was taken 
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Fig. 7. A comparison of the temperature structures of two red gi- 
ant atmospheric models. The top panel displays the temperatures 
as a function of logj()(Pgas). The solid line is a spherical model 
with the parameters L - 3690 Lq, M - 1 M© and R -\66Rq. 
These values were chosen to match the parameters of the plane- 
parallel model having Teff - 3500 K and \ogg - 0.0, computed 
using Atlas_ODF, shown by the dashed line. The bottom panel 
shows the percentage difference between the temperatures of the 
two models, also as a function of the gas pressure. 



Fig. 8. A comparison of the structures of two spherical red gi- 
ant atmospheric models, both equivalent to T^s - 3500 K and 
log^ - 0.0. The top panel displays the temperatures as a func- 
tion of logioCPgas)- The solid line is a spherical model having 
L = 3690 Lq, M = 1 M0 and R = 166 Rq. The dashed 
line, which is nearly coincident with the solid line, represents 
the model having the parameters L - 2952 Lq, M - 0.8 Mq 
and R - 148 Rq. The bottom panel shows the percentage dif- 
ference between the temperatures of the two models, also as a 
function of the gas pressure. The sense of the differences is the 
less massive and luminous model minus the more massive and 
luminous model. 



from the grid on Kurucz's web site, that model served only as 
the starting point for computing the model structure using our 
Atlas_ODF code to ensure that it reflects the same numerical 
routines. In particular, we used the Rybicki routine for the ra- 
diative transfer of both the plane-parallel and the spherical mod- 
els so that the resulting differences must come from the atmo- 
sphere's geometry, and not the method of solution. 

The spherical model is cooler than the plane-parallel model 
throughout most of the atmosphere, and it becomes progres- 
sively cooler with increasing height. The distance from to 
r(TR - 2/3) - Rt is 2.12 x 10^ km, giving an atmospheric exten- 
sion, defined in section 3.1, of 0.18. In the plane-parallel model 
the corresponding distance d(l) - d{TR = 2/3) - 1.85 x 10^ km. 
Deep in the atmosphere the spherical model becomes systemat- 
ically hotter than the planer-parallel model as the core makes a 
greater contribution. 

Because of the degeneracy of the atmospheric parameters, 
we tried another combination of luminosity, mass and radius 
that also match T^ff - 3500 K and \ogg = 0.0, namely, L, - 
2952 Lq,M = 0.8 Mq and R = 148 Rq. The comparison of 
the two spherical models, both computed with SAtlas_ODF, is 
shown in Fig. [8] The structures are so similar that they seem 
identical in the top panel. In the bottom panel, where the dif- 
ferences in the structures are plotted, it is easier to see that the 
less massive and luminous star has a slightly steeper temperature 
profile, as is expected. 



5. Comparison with other programs 

The Phoenix program ( Hauschildt et al.lll999h can also compute 
LTE, line-blanketed, spherically extended models, and a com- 
parisons with those models is appropriate. The Phoenix web site 
(http://www.hs.uni-hamburg.de/EN/ForfThA/phoenix/index.html) 
contains the NG-giant grids, in which the model that is closest 
to the examples used above is the one with Teff = 3600 
K, \ogg - 0.0 and M - 2.5 Mq. To compare with 
this model, we have computed a spherical model with 
L = 10324 Lq, M = 2.5 Mq and R = 262 Rq, again starting 
from the same plane-parallel model with Teff = 3500 K and 
\ogg - 0.0 that we used earlier. The comparison is shown in 
Fig. |9] Now the diff'erences are somewhat larger than in the 
previous comparisons, which is to be expected because the 
detailed calculations are nearly totally independent. Overall, 
however, the comparison is very similar, showing that the two 
models have essentially the same structures, although we note 
that the NextGen model has a temperature bulge compared with 
our model in the pressure range -1 < logiyPgas < +1.5. We 
have not observed this kind of feature in the comparisons we 
have made with the Kurucz models or between our spherical 
and plane-parallel models. 

A note about the relative run times of the two codes: the time 
per iteration running S Atlas_ODF on our single-processor desk- 
top workstation is just 5% of the time per iteration given in the 
header files of the NG-giant model. 
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Fig. 9. A comparison of the temperature structures of an 
SAtlas_ODF atmosphere (soHd line) with the structure of a 
model from the NG-giant grid computed with the Phoenix code 
(dashed line). The SAtlas model has the atmospheric parame- 
ters L = 10324 Lq, M = 2.5 M© and R = 262 Rq, correspond- 
ing to the NextGen parameters T^s - 3600 K, \ogg - 0.0 and 
M = 2.5 Mq. 

The MARCS program jGustafsson et al.1 l2008h is another 
well established code that has the ability to compute LTE, line- 
blanketed, spherical model atmospheres. From the MARCS web 
site {http://marcs.astro.uu.se/) the model with parameters T^s - 
4000 K, log^ = 0.0 and M - 1 Mq is closest the the ex- 
amples we have been using. This model also has solar abun- 
dances and a microturbulent velocity of 2 km/s. The header lines 
in the model gives the spherical parameters L - 6390 Lq and 
R = 1.1550 X 10'^ cm = 166 Rq. MARCS defines the ra- 
dius at tr = 1.0, not at tr — 2/3 that we use, but this is a 
small difference. Therefore, we have started from the model with 
Teff = 4000 K, log^ - 0.0 and microturbulence = 2 km/s in 
the same grid (/grids/gridpQ0od£new/apQQk2od£new. dat) 
used earlier, and we have computed a spherical model with the 
luminosity, mass and radius of the MARCS model. The com- 
parison is shown in Fig [TO] The models agree very well in the 
range logigPgas < -1-1.5, where the NextGen model displayed 
a temperature bulge. But, the spherical MARCS model appears 
to have a pressure inversion at T > 6500 K, something that is 
not present in the comparison with the Phoenix model. Overall, 
however, the structures are in substantial agreement. 

6. Conclusions 

We have modified the robust, open-source, plane-parallel model 
atmosphere program Atlas to treat spherically extended geom- 
etry. The resulting spherical code, SAtlas, which is available 
in both opacity distribution function and opacity sampling ver- 
sions, was used to compute several test models. At high surface 
gravity the spherical model structure is essentially identical to 
the plane-parallel model structure. At low surface gravity, the 



Fig. 10. A comparison of the temperature structures of an 
SAtlas_ODF atmosphere (solid line) with a spherical MARCS 
model (dashed line). The SAtlas model has the atmospheric pa- 
rameters L - 6390 Lq, M = 1.0 Mq and R - 166 Rq, corre- 
sponding to the MARCS parameters Teff = 4000 K, log^ - 0.0 
andM = 1.0 M©. 

SAtlas models agree very well with the spherical model struc- 
tures computed by Phoenix and by MARCS. The SAtlas pro- 
gram, which runs easily on a desktop workstation, offers a viable 
alternative for modeling the atmospheres of low surface gravity 
stars. 

As an example of the utility of SAtlas, we have used it 
to compute more than 2500 models to create model cubes 
with fine parameter spacing covering the specific L», M» 
and /?, values needed for an analysis of the optical inter- 
ferometry of just three stars (Neilson & Lester submitted). 
We are in the process of computing more models for our 
own application. These codes and models are available at 
http : //www . astro . utoronto . ca/~lester/Prograins/. 
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